www.gusucode.com > 溷沌分析工具箱 - Chaotic Systems Toolbox > 混沌分析工具箱 - Chaotic Systems Toolbox\lorentz.m
function [x,y,z]=lorentz(n,level,s,r,b,x0,y0,z0,h) %Syntax: [x,y,z]=lorentz(n,level,s,r,b,x0,y0,z0,h) %_________________________________________________ % % Simulation of the Lorentz ODE. % dx/dt=s*(y-x) % dy/dt=r*x-y-xz % dz/dt=x*y-b*z; % % x, y, and z are the simulated time series. % n is the number of the simulated points. % level is the noise standard deviation divided by the standard deviation of the % noise-free time series. We assume Gaussian noise with zero mean. % s, r, b, and s are the parameters % x0 is the initial value for x. % y0 is the initial value for y. % z0 is the initial value for z. % h is the step size. % % % Notes: % The time is n*h. % The integration is obtained by the Euler's method. % % % Reference: % % Lorentz E N (1963): Deterministic nonperiodic flow. Journal of the % Atmosphairic Sciences 20: 130-141 % % % Alexandros Leontitsis % Department of Education % University of Ioannina % 45110 - Dourouti % Ioannina % Greece % % University e-mail: me00743@cc.uoi.gr % Lifetime e-mail: leoaleq@yahoo.com % Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421 % % 16 Nov 2001 if nargin<1 | isempty(n)==1 n=500; else % n must be scalar if sum(size(n))>2 error('n must be scalar.'); end % n must be positive if n<0 error('n must be positive.'); end % n must be an integer if round(n)-n~=0 error('n must be an integer.'); end end if nargin<2 | isempty(level)==1 level=0; else % level must be scalar if sum(size(level))>2 error('level must be scalar.'); end % level must be positive if level<0 error('level must be positive.'); end end if nargin<3 | isempty(s)==1 s=16; else % s must be scalar if sum(size(s))>2 error('s must be scalar.'); end end if nargin<4 | isempty(r)==1 r=45.92; else % r must be scalar if sum(size(r))>2 error('r must be scalar.'); end end if nargin<5 | isempty(b)==1 b=4; else % b must be scalar if sum(size(b))>2 error('b must be scalar.'); end end if nargin<6 | isempty(x0)==1 x0=0.1; else % x0 must be scalar if sum(size(x0))>2 error('x0 must be scalar.'); end end if nargin<7 | isempty(y0)==1 y0=0.1; else % y0 must be scalar if max(size(y0))>2 error('y0 must be scalar.'); end end if nargin<8 | isempty(z0)==1 z0=0.1; else % z0 must be scalar if max(size(z0))>2 error('z0 must be scalar.'); end end if nargin<9 | isempty(h)==1 h=0.01; else % h must be scalar if max(size(h))>2 error('h must be scalar.'); end % h must be positive if h<0 error('h must be positive.'); end end % Initialize y(1,:)=[x0 y0 z0]; % Simulate for i=2:n ydot(1)=s*(y(i-1,2)-y(i-1,1)); ydot(2)=r*y(i-1,1)-y(i-1,2)-y(i-1,1)*y(i-1,3); ydot(3)=y(i-1,1)*y(i-1,2)-b*y(i-1,3); y(i,:)=y(i-1,:)+h*ydot; end % Separate the solutions x=y(:,1); z=y(:,3); y=y(:,2); % Add normal white noise x=x+randn(n,1)*level*std(x); y=y+randn(n,1)*level*std(y); z=z+randn(n,1)*level*std(z);